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We develop a path integrals approach for analyzing stationary light propagation appropriate for 
photonic crystals. The hermitian form of the stationary Maxwell equations is transformed into a 
quantum mechanical problem of a spin 1 particle with spin-orbit coupling and position dependent 
mass. After appropriate ordering several path integral representations of a solution are constructed. 
One leaves the propagation of polarization degrees of freedom in an operator form integrated over 
paths in coordinate space. The use of spin 1 coherent states allows to represent this part as a 
path integral over such states. Finally a path integral in transversal momentum space with explicit 
transversality enforced at every time slice is also given. As an example the geometrical optics limit 
is discussed and the ray equation is recovered together with the Rytov rotation of the polarization 
vector. 
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Introduction. Description of stationary light fields in 
photonic crystals can be formulated, cf. as solu- 

tions of an eigenvalue equation. 



V X 



— -VxH= - 
e(x) V c 



H 



(1) 



where H is the magnetic field and e(x) is the space de- 
pendent dielectric constant. This equation is sometimes 
referred to in the literature as the master equation. The 
operator = V x -t^-Vx is hermitian and positive defi- 
nite. This formulation allows direct applications to pho- 
tonic crystals of many powerful techniques developed in 
quantum mechanics. 

Our goal is to follow on this development and formu- 
late the path integral representation of solutions of this 
equation. Within the scalar approximation 
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(2) 



works in this spirit have already been presented in the 
past [3, 0, H, S, l3l , but here we deal with the full vector 
version. To the best of our knowledge this have not yet 
been done. 

Let us consider an auxiliary time dependent 
Schrodinger like equation 



idrC{x,T) = eC(x,T) 



(3) 



with fictitious time parameter r. This equation has the 
formal solution 



C(x,r) = exp -ire C(x,0) 



(4) 
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from which we can recover a solution of the original mas- 
ter equation by using 



H(x) 



lim 

?)^0 



dT 



(5) 



Note that for any uj ^ the solutions H(x) of Eq. ([T]) are 
transversal, V • H = 0. The solutions C (x, r) of ^ are 
not. However since 9r [V • C(x, r)] = we have that the 
above transversality condition is satisfied also by H(x) 
given by ([5|). 

The Dirac notations. The Hilbert space of all complex 
vector functions H(x) is equivalent to the Hilbert space 
of a spin 1 particle in quantum mechanics. The spin 
operators are hidden in the vector product signs which 
we make explicit by using the antisymmetric tensor (and 
the summation over repeated indices convention) 
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e(x) 



{-idi)Hm., 
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The operator G) is therefore 

6= (p-S)z;(x)(p.S) 

with p — — iV , 

v{x) ~ l/e(x) 



and S - the spin 1 matrix vector with components 

['5'j]ifc = ^i^ikj ■ 

It is convenient to introduce the Dirac notations 



e\H) 

{x,i\Q\H) 
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e(x) 



i 1,2,3 
V X H (x) 
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In this notation equation ((4]) is 



(x,z|C(r))=^ / dx'(x,^|e— «|x',j)(x',j|C(0)) (9) 

Ordering the " Hamiltonian" . Wc intend to write 
the path integral expression for the propagator 
(x, z|e~*'^®|x', j). FoUowing the usual procedure Q we 
rewrite O in the ordered form placing the p operators to 
the right of x's, 



e = -z(Vu • S)(p • S) + v{x)p'^'Pt 
with Vt ^ the transversal projection operator 



(10) 



(11) 
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We are interested only in how the transverse functions 
C(x, r) propagate for which Vt\C) = |C). Accordingly 
we can drop Vt in Eq- (flOll and define a reduced operator 



(12) 



efl = ^;(x)p2-i(Vi;-S)(p.S) 



where we used (SiSj)r 



Note that the operators ([7|) and are equivalent when 
acting on transversal functions. In general one can show 
that Q ~ QrVt- Note that is hermitian in the trans- 
verse subspace. 

Path Integral - The Operator Version. We are now in a 
position to use the standard time slicing process to con- 
struct the path integral. In this we first choose to keep 
the vector part of the evolution at the operator level in- 
serting complete (coordinate and momentum) states only 
for the spatial part. As a result we obtain a functional 
integral over the positions and momenta of spatial paths 
with the integrand containing the evolution operator of 
the vectorial part of the field for each path 



-*^®|x') = lim 

At^O 



dxidpi . . . dp 



N 



N 



3N 



(2^) 

^[x(T),p(T)]exp 



-Q gzp,.(x,-x,_i) _ -^^^^ p2j ^ _ Ar(Vwj • S)(pj • S) 



(p • 9i-x — tip ) dr 



r<^exp[- / (VvS)(p-S)dT] 



(13) 



where T stands for time ordering, Ar = t/N, x(0) = x', 
x(r) = X. This functional integral has position depen- 
dent mass m(x) = l/2u(x) = e(x)/2 and spin orbit cou- 
pling terms which enter the time ordered exponential. 

The Spin Coherent States Version. The operator part 
of the above path integral, i.e. the time ordered expo- 
nential can be further developed using the spin coherent 
states. These states are generated, cf. [13, [iM [H, [l7| . 
by rotations of one of the eigenstates of Sz, i.e. \l,fJ,), 
/i = 1,0, — 1. In the cartesian basis used above these 
states are 




1,0) = 




1,-1) = 




The most common spin coherent states are given by 



\n) = |0,0) = 



"11,1) 



(14) 



They are eigenstates of the spin component in the appro- 
priate direction, ft ■ S|0) = jSl), and form an over com- 
plete set which has a useful " resolution of unity" property 



47r 



dn\n){n\ , dn = siii9ded(j). (15) 



Going again through the time slicing process and insert- 
ing (fT5|) between the slices one finds the propagator be- 
tween two spin coherent states 



x/,^2/|e-'^««|x„r!,) = / 7^[x,p] / i^Pexp 



i I (p-drx + cos0d^(j)-vp^ + i{n\{Vv-S){p-S)\n))dT 



(16) 



All that is left is to evaluate the matrix element Using SiSj — ^{[Si, Sj]+{Si, Sj}) together with the com- 

{n\ivv ■ s)(p • s)\n) = Y,i'^^v)pJ{^\s^s,\n) 
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mutation relations of spin operators, the identity 
valid for spin 1 matrices, [3| as well as 



(17) 



cf. [17|, we obtain 



(18) 



and finally 



{n\{Vv ■ s)(p • s)\n) = -n ■ (Vw x p) 



(19) 



Explicit Transversal Projection. The above path inte- 
gral expressions for the exact propagator will propagate 
only transverse vector functions despite unrestricted in- 
tegrations at every time slice. It may be desirable how- 
ever, especially in making approximations, to have path 
integral expressions with explicit transversal projection 
enforced at every time step. This can be achieved by in- 
serting the projection operator (jlip at every time slice. 
Moreover one must start with a transverse state, the sim- 
plest of which is A exp(ip • x) with A • p = 0. It is natural 
then to take also the final state to be a transversal plane 
wave. Using the spin coherent states this means working 
with (py, riy|e~''^®"|pi, r^i) and imposing the transver- 
sality conditions 



n 



± 



Pi,/ 

\Pi,f\ 



(20) 



under which VT\Pi.f,^i,f) = \PiJ:^i,f) ■ Once one uses 
such transversal initial and final states and due to the 
fact that the transversality is conserved by the infinitesi- 
mal transformations in between infinitesimal time slices, 
the projection operator insertion will not alter the prop- 
agation. Expression (jl3[) will then become (in the limit 
At 0) 



(p/|e— Q«|p,> = 



J -D [p, x] cxp[— i J (x • drP + wp^) dr] 
n[l-Ar(vw, -S) (p, -S'' ''^^ 



The Saddle Point Approximation - The Geometrical 
Optics Limit. Geometrical optics is a short wavelength 
expansion. The vacuum wavelength of light must be 
short with respect to the scale w/|Vu| over which the 
dielectric function is changing. Accordingly we treat the 



spin dependent part of Eqs. p3|) and p6|) as slowly vary- 
ing and evaluate it on the saddle point of the functional 
Sq = Jq {p ■ drX — vp'^)dT. The corresponding Euler- 
Lagrange equations are 



drX = 2vp ; drP = — p Vw 



(21) 



The conserved " energy" is E = vp^ which we will use in 
order to reparametrize 

dr^2vVEdt (22) 

and obtain the standard differential equation of the ray 

d^x „fn^\ 

Once the geometrical ray (or rays) satisfying the appro- 
priate boundary conditions (e.g. x(0) — x',x(r) — x) is 
found the propagator along each ray is given 



K = T- 



-/o-(V«,-S)(p,.s)d 



T|giSo(x,x') j 



D[^,p\ 



(24) 

where the subscript r means evaluation at the ray val- 
ues which also applies for S'o and (5^ 5*0. The Gaus- 
sian integral over expi(5^S'o/2 is evaluated in a stan- 
dard way 0, giving the Van-Vleck determinant factor 
^(l/27r)3det{9^^,5o(x,x')}. The evolution of the vec- 
tor degrees of freedom along a geometrical ray is given 
by the equation 



dr 



V«,-S Pr-S |A) 



(25) 



as follows from the time ordered operator in ((24l) . In the 
usual vector representation this is just an evolution of a 
vector A 



d\ 

d^ 



VVr X {pr X A) 



(26) 



with time dependent Vu^ = Vw(xr(T)) and Pr(T) deter- 
mined by the ray trajectory Xr(T), Pr(T). We will now 
show that the orientation of A relative to the ray is gov- 
erned by the Rytov equation for the polarization while 
its magnitude fits one of the components of the energy 
flow equation along the ray, cf. 3, 20| . 

It is convenient to switch from "time" r to the length 
parameter s and write (|26p as 



dX 

ds 



Vn 



X (T(,s) X A) 



(27) 



where T(s) = dx/ds is the tangent to the curve. Consid- 
ering also the normal and the binormal vectors 



N{s) 



dT{s) 
ds 



dT{s) 



ds 



B(s) = T(s) X N(s) 
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one can expand 

A = AtT + AatN + AbB. (28) 
We transform the equation ([23]) of the ray from r to s 

(29) 



-T- n— I = Vn 

as \ as 



and use the Frenet equations, describing the differential 
geometry of curves, [23| 



(30) 



where k{s) = \dT/ds\ is the curvature and r (s) = B • 
d'N/ds is the torsion (not to be confused with the time 
parameter r). We can rewrite Eq. (P7|) as 









(T 


-K 




N 










d 

ds 







it) 










(31) 



The At (s) component is conserved along the curve and 
we will take it to be zero (a transversal solution). We 
will rewrite the remaining two equations as equations for 
(fi ~ the angle between A and N and A^ = A^ + A^- the 
square length of A. We obtain 



dip 
ds 



dA2 
"d7 



dln[r 



ds 



(32) 



The first equation is the Rytov equation from geomet- 
rical optics governing the rotation of polarization along 
an optical ray, cf.. Fig. [1] The solution of the second 
equation is simply X = n{s) Xq. 
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FIG. 1: Rytov rotation of the field relative to the ray. T, N 
and B are respectively tangent, normal and binormal vectors. 
H is the magnetic field vector, x is the geometrical ray. 

In the standard derivations of the geometrical optics 
from the Maxwell equations, cf., 0,11, 

one finds the 



following equation for the change with respect to the 
length parameter of the field amplitude square along a 
geometrical ray 

'^Ink^l^lAI^ (33) 



d 
ds 



ds 



The two terms in the right hand side are clearly asso- 
ciated with the two contributions which we derived in 
the geometrical optics approximation to the amplitude 
of the propagator. The first term, which involves the 
second spatial derivatives of the action, is related to the 
Van-Vleck determinant contained in the last term in Eq. 
([24]) . The second term is identical with what we obtained 
for A2 in Eq. ([32]). 
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